function Replicate(DATA, PARAMS, OUT)

    close all
    format long g
    
    tolerance = 1e-12;
    
    addpath('helpers')
    addpath('rcnl_demand')
    addpath('rcnl_supply')
    
    
    [df, dfull, cdindex] = load_data_full(DATA);
    [xlinear, ~, xnonlin, fe, ~, zdemand, ~, ~, xid] = final_variables(df, 'og', 'og');
    
    s_jt     = df.s_jt;
    p_jt     = df.p_jt;
    prodid   = df.prodid;
    market   = df.cdid;
    cluster  = df.cityid;
    dresults = struct;
    Pi       = [0.0002; 0.035; 0.004]; % OG starting values from MW
    rho      = 0.5;
    par_rest = 1:3;
    varargin = {};

    coefs_sim = struct;
    coefs_sim.dist_mis = nan;

    try
        [dresults, IVDtilda] = estimate_rcnl( ...
            Pi,      ...
            rho,     ...
            s_jt,    ...
            p_jt,    ...
            prodid,  ...
            xlinear, ...
            xnonlin, ...
            zdemand, ...
            zdemand, ...
            fe,      ...
            market,  ...
            dfull,   ...
            cluster, ...
            "RCNL",  ...
            par_rest, ...
            coefs_sim, ...
            tolerance, ...
            varargin{:});
        dflags = 0
    catch err
        fprintf('%s: %s\n', err.identifier, err.message);
        dflags = 1
        IVDtilda = NaN;
    end
    
    if dflags == 0
        sresults = estimate_rcnl_supply( ...
            df,                          ...
            cdindex,                     ...
            dresults,                    ...
            dresults.derMat,             ...
            IVDtilda,                    ...
            'og'                         ...
        );
    else
        sresults = struct;
    end
    
    ddedup = dfull(cdindex, :);
    
    outfile = fullfile(OUT, 'Replicate.mat');
    
    save(outfile,   ...
        'df',       ...
        'ddedup',   ...
        'cdindex',  ...
        'dresults', ...
        'dflags',   ...
        'sresults'  ...
    );
    
    assertion = max([abs(dresults.rho - 0.8299),  ...
                abs(dresults.theta(1) + 0.0887), ...
                abs(dresults.Pi(1) - 0.0007),    ...
                abs(dresults.Pi(2) - 0.0143),    ...
                abs(dresults.Pi(3) - 0.0043)]);
    
    if assertion < 0.0001
        disp(' ')
        disp('Parameters estimated correctly.');
    else
        disp(' ')
        disp('Parameters depart from the ones reported in the paper.');
        disp(['rho = ',num2str(dresults.rho)]);
        disp(['alpha = ',num2str(dresults.theta(1))]);
        disp(['Pi_1 = ',num2str(dresults.Pi(1))]);
        disp(['Pi_2 = ',num2str(dresults.Pi(2))]);
        disp(['Pi_3 = ',num2str(dresults.Pi(3))]);
    end
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Compute auto-filled scalars %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    addpath('../../analysis_cluster/code/helpers')
    addpath('../../analysis_cluster/code/Sim')
    
    [df, dfull] = load_data_full(DATA);
    df_covsize = df;
    [~, ~, coefs]  = load_parameters(fullfile(PARAMS, 'Replicate.mat'));
    
    % Scalars: R2 of calendar month fixed effect on the seasonal month.
    
    sdate    = f_selectfe(df.dateid, df.dateid, 1);
    sdatem   = accumarray(df.dateid,df.montid,[],@mean); 
    dateD_fe = [0;coefs.sigma_D(sdate+max(df.prodid)-1)];
    montD_fe = accumarray(sdatem(2:end),dateD_fe(2:end),[],@mean);
    montD_fe = montD_fe(sdatem);
    
    mdl = fitlm(montD_fe,dateD_fe);
    
    R2 = 100 * mdl.Rsquared.Ordinary;
    
    myList = struct;
    myList.RSquared = round(R2, 1);
    save_scalars(myList, '../output/product/scalars/R2', '', false);
    
    % Scalars: ratio of the variance of the simulated income to the variance of MW's specification and
    % the ratio of the variance of the simulated income to the variance of the market-level income.
    
    C             = size(dfull,2); 
    cityyearid    = grp2idx(df.yearid*100+df.cityid);
    mean_inc      = accumarray(cityyearid,df.minc,[],@mean);
    income_draws  = cell2mat(arrayfun(@(x) accumarray(cityyearid,dfull(:,x),[],@(x)mean(x)),1:C,'un',0))+mean(mean_inc);
    income_norm   = income_draws./mean_inc;
    income_norm   = income_norm(:);
    income_sample = datasample(income_norm,C)';
    minc_new      = coarsen_income(mean_inc);
    df.minc       = minc_new(cityyearid);
    
    dfull_sim     = repmat(income_sample,size(df.id2)).*df.minc;
    dfull_sim     = dfull_sim - mean(mean(dfull_sim));
    
    Var_psi     = accumarray(df.cdid,var(dfull,[],2),[],@mean);
    Var_psi_sim = accumarray(df.cdid,var(dfull_sim,[],2),[],@mean);
    
    Variance_ratio = 100 * mean(Var_psi_sim)/mean(Var_psi);
    
    Variance_ratio_market = 100 * var(accumarray(df.cdid,mean(dfull_sim,2),[],@mean))/var(accumarray(df.cdid,mean(dfull,2),[],@mean));
    
    myList = struct;
    myList.VarianceRatio = round(Variance_ratio, 1);
    myList.VarianceRatioMarket = round(Variance_ratio_market, 1);
    save_scalars(myList, '../output/product/scalars/Variance_ratio', '', false);


    df_covsize.inc_mean = coarsen_income(df_covsize.minc);
    df_covsize.prodid_string = string(df_covsize.prodid);
    [groupID, groupNames] = findgroups(df_covsize.cdid);
    prod_list = strings(height(df_covsize), 1);
    for i = 1:max(groupID)
        prod_ids = join(df_covsize.prodid_string(groupID == i), ', ');
        prod_list(groupID == i) = prod_ids;
    end
    df_covsize.prod_list = prod_list;
    df_covsize.X_id = strcat(string(df_covsize.montid), '_', string(df_covsize.prod_list), '_', string(df_covsize.inc_mean));
    [unique_categories, ~, category_idx] = unique(df_covsize.X_id, 'stable'); 
    df_covsize.X_id = category_idx;

    unique_comb = unique(df_covsize(:, {'cdid', 'X_id'}));

    writetable(unique_comb, '../output/unique_combinations.csv');
    
    end

